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Gravitational waveforms and fluxes from extreme mass-ratio inspirals can be computed using 
time-domain methods with accuracy that is fast approaching that of frequency-domain methods. 
We study in detail the computational efficiency of these methods for equatorial orbits of fast spin- 
ning Kerr black holes, and find the number of modes needed in either method — as functions of the 
orbital parameters — in order to achieve a desired accuracy level. We then estimate the total com- 
putation time and argue that for high eccentricity orbits the time-domain approach is more efficient 
^v^j \ computationally. We suggest that in practice low-m modes are computed using the frequency- 

domain approach, and high-m modes are computed using the time-domain approach, where m is 
the azimuthal mode number. 
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Compact stellar mass objects emit gravitational waves in the good sensitivity band of the planned Laser Intcrfcr- 
^Jy ometer Space Antenna (LISA) during their last year of inspiral into a supermassive black hole. The gravitational 
waveforms of EMRIs, extreme mass ratio insiprals, are therefore of interest, and when observed may shed light on 
the spacetime of the central black hole, including testing the Kerr hypothesis in addition to being a sensitive tool 
to measure the central black hole's parameters. This possibility makes such sources of gravitational waves extremely 
interesting, despite the relatively low number of sources expected during the lifetime of the LISA mission. Specifically, 
tens to hundreds of such sources are expected to redshifts of z ~ 0.5-1 0. In addition, the solution of the two body 
problem in general relativity remains challenging, and EMRIs orbits and waveforms correspond to its solution in the 
extreme mass ratio limit. 
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PACS numbers: 04.25.Nx, 04.30.Db, 04.30.- 



I. INTRODUCTION AND SUMMARY 



Construction of theoretical templates is important both for detection of EMRIs gravitational waves and for accurate 
parameter estimation. Numerical waveforms can be constructed using the frequency-domain (FD) or the time- 
domain (TD) approaches. The former approach has been developed to very high accuracy, and is considered robust 
and accurate. On the other hand, advances to the TD approach have been hindered first by the success of the FD 
approach d, 0| , and by the crudity of the initial attempts to evolve numerically the fields coupled to a point-like 
source with the Teukolsky equation [f| Q . The breakthrough in the accuracy of TD solutions of the inhomogeneous 
Teukolsky equation, i.e., the 2+1D solution of the Teukolsky equation coupled to a point mass, was recently achieved 
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For the first time, it was shown that TD calculations can be as accurate as FD calculations. The TD method 
jj was improved with the introduction of the "discrete delta" model of the source Q and an appropriate low pass 
filter that makes the discrete delta useful also for eccentric or inclined orbits @ . Specifically, correlation integrals of 
gravitational waveforms done for the same system in the FD and TD approaches show that the two agree to a high 
level 0. One may therefore argue that the two methods are comparable in the results they are capable of producing. 
We therefore contend that the viewpoint that the TD solution of the inhomogeneous Teukolsky equations is far from 
being competitive from the FD solution can no longer be supported. However, we believe that one should not seek 
competition of the two approaches, but rather how they complement each other, as either method has non-overlapping 
strengths. In order to achieve this goal one needs to compare the computational efficiency of the two approaches. 

The question of the efficiency and the computation time with which the results are obtained remains an open question 
though. The common wisdom is that the FD computation is more effective computationally than its TD counterpart: 
FD approaches are particularly convenient when the system — and the emitted gravitational waves — exhibits a 
discrete set of frequencies. Indeed, as shown by Schmidt [l(| and by Drasco and Hughes [TT| . all bound Kerr orbits 
have a simple, discrete spectrum of orbital frequencies. However, generic orbits and in particular high-eccentricity 
orbits, although in principle amenable to a Fourier decomposition and a FD construction of the waveforms, require 
the summation of many terms in the Fourier series. This problem limits the accuracy and increases the computation 
time in FD calculations. While this statement, or similar ones, appear in the literature [3j, it has not been quantified. 
The motivation of this paper is to study the relative computational efficiency of FD and TD codes for the solution 
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of the inhomogeneous Teukolsky equation. Specifically, we study the question of how much computational time is 
needed to find the fluxes of energy and angular momentum to infinity from an eccentric and equatorial orbit around a 
fast spinning Kerr black hole. We restrict the analysis here to equatorial orbits. Analysis of equatorial orbits is non- 
trivial, and teaches much of the method and properties also of non-equatorial orbits. Analysis of non equatorial orbits 
increases considerably the volume of the parameter space, and we leave its study to the future. We do, nevertheless, 
comment on the analysis of generic orbits. 

The estimation of the total needed computation time for TD codes is relatively simple: after the desired accuracy 
level is set, one needs to estimate the number of azimuthal m modes (associated with the <f> orbital angular frequency) 
required for the total sum over m modes to achieve the desired accuracy, and then compute each m mode at the 
same accuracy level. The number of needed m modes is not hard to estimate, because the partial fluxes approach 
a geometric progression in m for large m values, with the asymptotic factor between two successive partial fluxes 
depending on the eccentricity of the orbit. After the needed number of m modes is determined, one can use the TD 
code of @ (possibly with the improvements included in to calculate the partial fluxes and sum over them. Using 

the asymptotic geometric progression structure one may also estimate the error associated with the summation over 
the partial fluxes, and verify that the desired accuracy level is indeed achieved. This is done in Section [Til The TD 
calculation of the individual partial fluxes can be done more efficiently than a straightforward calculation of the fields 
to very late times and great distances, that are required for an accurate estimate of the fluxes at infinity. Specifically, 
one can make use of the peeling properties of the Weyl scalars at great distances, and fit the fields at finite distances 
to that behavior. One may then use the fitted behavior to extract the behavior at infinity. This is done in Appendix 
|A1 We find this method to be cleaner and better motivated physically than the fit done in [8|, that uses a general 
two-parameter fit function of a form which is not strongly motivated physically. 

The estimation of the total needed computation time for FD codes starts similarly to the TD analysis: specifically, 
one needs to first determine the number of m modes needed for the determined accuracy level. This can be done 
in a similar way to how it is done in the TD. The calculation of each m mode in the FD approach requires the 
summation over a number of k,n,£ modes, that correspond to the radial and angular (about the equatorial plane) 
orbital angular frequencies and the radiative multipole t. Because we focus attention on equatorial plane orbits, the 
n,£ modes trivialize, and in practice only the k modes are of importance. We find in Section lllll that the number of 
k modes that one needs to sum over to achieve the desired accuracy level increases with the eccentricity and with the 
corresponding value of m in an intricate manner. 

Next, in Section [TVl we estimate the computation time of each k mode, and find it is a function of k. The behavior 
of the computation time of the k mode is found to have a rather intricate dependence of the value of k. We use this 
behavior to approximate the computation time over a sum of k modes up to some fc max , and argue it is approximately 
quadratic in fc max for high fc max values. Then, in Section [V] we find the total computation time for all the FD modes 
on the same machine on which we perform the TD computation, and compare the two. We find that the FD code 
is more efficient at low m values (for a twofold reason: because those require few k modes, and for low m values 
each k mode takes less time to compute), but the computation time increases rapidly with the value of m. We find 
the growth rate of the computation time with the mode number m to increase as a power of m, with the value of 
the exponent increasing with the eccentricity. Finally, we estimate that for generic orbits the two methods become 
comparable already for moderately high values of the eccentricity, in the range of e ~ 0.6-0.7. Higher eccentricity 
orbits are more efficiently calculated using the TD approach. Even when the calculation of the sum of all m modes 
is done more efficiently using one method over the other, one may still compute the low m modes using the FD 
approach and the high m modes using the TD approach. Such a hybrid method may prove to be the most efficient 
computationally. One should also consider the fact that the total needed number of k modes is an empirically found 
number. When the full parameter space is mapped, one may tabulate the numbers of required modes as functions 
of the orbital parameters. Until this is done, extra computation time is needed to find the number of needed modes, 
including the computation of some modes that are found a posteriori to be unneeded. No similar problem occurs in 
the TD approach. 



II. SUMMING OVER THE rn MODES 



To obtain the sum over all m modes we need to obtain results for high m-numbers, and find a way to determine (i) 
how many m modes we need to sum over to get the total flux to a certain pre-determined accuracy, and (ii) estimate 
the error in neglecting all the higher m modes. In the case of circular and equatorial orbits, Finn & Thorne show 
that (see [12j for more details and for definitions) 

■ _ 2(m + l)(m + 2)(2m + 1)! m 2m + 1 „ 2ft2+2m/3 d 
m ~ (m-l)[2"»m!(2m + l)H]2 V toom 
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which has the nice property that E m+ i/E m — ► const as to — > oo. (Notably, this property depends on the factor of 
may therefore approximate the sum over infinitely many modes by taking the series to be a geometric 

progression. 

For eccentric orbits we no longer have the Finn-Thorne formula, but we can do numerical experiments to test 
whether the same results holds also for such orbits. We first show in Table Q] the average fluxes obtain for a number 
of eccentric orbits (with a fixed value of the semilatus rectum p) as a function of the mode m. The data presented 
is Table U is intentionally coarse, and is not intended to be more accurate than at the 5-10% level. Below, we also 
present similar data with higher accuracy, where the latter are needed. We plot in Fig. [T] E m as a function of m 
for several eccentric orbits. We see that for all values of e, the drop off rate is exponential in to for large values of 
to. Notably, the highest flux is in the fundamental mode, to = 2, and is highest for e ~ 0.5. This result is indeed 
expected: The radiated power in gravitational waves averaged over one period for a point mass \x in a Keplerian orbit 
around a Schwarzschild black hole of mass M is given by the quadrupole formula to be [l3| 

32 M 2 M 2 (M + M ) A 73 2 37 4 
[ ' 5 a 5 (l-e 2 ) 7 / 2 ^ + 24 96 6 

Recalling that the semimajor axis a is related to the semilatus rectum p by p = a(l — e 2 ), we find that for fixed p, 

(P) ~ (1 - e 2 ) 3 / 2 ( 1 + -e 2 + —e 4 ) , 

which has a maximum at e = 0.465 fl4| . As the partial flux in any other mode is suppressed by over an order 
of magnitude compared with the fundamental mode (m = 2), this result for the total flux is carried over to the 
fundamental mode. Notably, the larger to, we find numerically that the larger the eccentricity for which most flux is 
obtained. 



TABLE I: Average (over one period) fluxes per unit mass extracted at r = 100M for a central black hole with a/M — 0.9, and 
a prograde eccentric orbit of semilatus rectum p = 4.64M, for various values of the eccentricity e. For each mode m, we show in 
bold print the mode that corresponds to the eccentricity for which the highest flux is achieved. Notice that the data presented 
in this Table are rather coarse, and are only intended to demonstrate the overall behavior. The accuracy level of data in this 
Table is at the 5-10% level. Fine details are studied below. 



M 


e = 0.1 


0.5 


0.7 


0.8 


0.9 


i 


1.03 x 10 -6 


1.32 x 10" 6 


1.19 x 10~ 6 


1.07 x 10 -6 


1.05 x 10" 6 


2 


6.96 x 10~ 4 


8.73 x 10 4 


7.84 x 10~ 4 


6.49 x 10~ 4 


3.32 x 10~ 4 


3 


1.50 x 10~ 4 


2.57 x 10~ 4 


2.68 x 10 4 


2.37 x 10~ 4 


1.30 x 10~ 4 


4 


4.12 x 10~ 5 


9.34 x 10~ 5 


1.10 x 10 4 


1.03 x 10~ 4 


5.86 x 10" 5 


5 


1.16 x 10~ 5 


3.59 x 10~ 5 


4.92 x 10 5 


4.61 x 10~ 5 


2.86 x 10~ 5 


6 


3.58 x 10~ 6 


1.50 x 10~ 5 


2.31 x 10 5 


2.28 x 10~ 5 


1.48 x 10~ 5 


7 


1.06 x 10~ 6 


6.48 x 10~ 6 


1.12 x 10~ 5 


1.16 x 10 5 


7.80 x 10~ 6 


8 


3.25 x 10~ 7 


2.83 x 10" 6 


5.56 x 10" 6 


6.07 x 10 6 


4.28 x 10" 6 


9 




1.27 x 10~ 6 


2.79 x 10~ 6 


3.21 x 10 6 


2.37 x 10~ 6 


10 




5.72 x 10" 7 


1.42 x 10~ 6 


1.72 x 10~ 6 


1.33 x 10" 6 


11 




2.65 x 10~ 7 


7.30 x 10~ 7 


9.31 x 10 7 


7.53 x 10~ 7 


12 




1.25 x 10~ 7 


3.63 x 10~ 7 


5.05 x 10 7 


4.19 x 10~ 7 


13 




5.50 x 10~ 8 


1.94 x 10~ 7 


2.76 x 10 7 


2.39 x 10" 7 



We next consider the ratio of the relative contribution of the m mode to the partial-sum (up to to) flux compared 
with the preceding mode, as a function of to. An asymptotic drop off of the flux corresponds to this ratio approaching 
a constant value as to — > oo. In Fig. [2] we plot this ratio as a function of to for various values of the eccentricity. 

Lastly, we extrapolate the curves in Fig.[5]to find the asymptotic value for the ratio, and plot it in Fig.[3]as a function 
of the eccentricity e. The best fit curve is given by R = 0.24 + 0.44 e with a correlation coefficient R 2 = 0.980438. 

As E m behaves asymptotically like a geometric progression, we can approximately sum over all modes (provided 
sufficiently many modes are calculated, so that the sequence of partial fluxes already converges approximately to the 
asymptotic behavior). Specifically, calculating the sequence E\, Ei, E3, ■ ■ ■ , E n -i, E n , we can calculate the partial 
sums S n := J2m=i Em, and the remainder is approximated by R n ~ E„R/(1 — R), where R is the asymptotic ratio 
evaluated above. One can then approximate the total flux by £Totai ~ S n + R n , and use A n — Rn/ -^Totai &s ci measure 
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FIG. 1: The flux in the m mode as a function of m for various values of the eccentricity e. The black hole's spin is a/M = 0.9 
and the orbit's semilatus rectum is p/M = 4.64. 
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FIG. 2: The ratio of the relative contribution of the m mode to the partial-sum (up to m) flux compared with the preceding 
mode, as a function of m, for various values of the eccentricity e. Same parameters as in Fig. [T] 



for the error. We can therefore answer questions such as "how many modes m do we need to sum over to obtain a 
desired accuracy level?" (assuming each mode individually has the desired accuracy level). We show this in Fig. |4j 
that shows contour curves at a fixed accuracy level on the number of m-modes-eccentricity plane. 

Based on Fig. [4] to determine the energy flux to 10% accuracy requires 4 modes at all values of the eccentricity 
e. We may therefore find the approximate total waveform by summing over the 4 modes of greatest energy flux. For 
a/M = 0.9, p/M = 4.64, and e = 0.1, we show in Figs. [5] and O the waveform obtained when summing over these 4 
modes, specifically m — 2, 3, 4, 5. Indeed, the phase information is captured well by the sum of these 4 modes, and is 
no longer changing significantly by addition of more modes. 
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FIG. 3: The asymptotic ratio R as a function of the eccentricity e. The individual error bars are shown, in addition to the best 
fit curve (solid) and 3<r confidence curves (dotted). 



The behavior of the m-modes — that corresponds to the angular frequency Cl<f, — is common to both time-domain 
and frequency-domain approaches. In the TD case, however, the study of the behavior of the m-modes summarizes all 
the modes that need to be considered to get the full waveform and the total fluxes radiated. In the FD case, however, 
one needs to consider also the modes that arise from the three different frequencies of the problem: the fc-modcs, 
corresponding to the radial angular frequency f2 r , and the n-modes, corresponding to the inclination angle oscillations 
with angular frequency r2#. (We note that various authors exchange the notation of the latter two frequencies, k «-» n.) 
Here, we concentrate on the equatorial plane, and defer discussion on motion outside the equatorial plane to a sequel. 

The FD code that we use is based on the code presented in Q, where more details on the code can be found in 
addition to descriptions of the tests done to check the code and the sources for errors. The FD numerical method is 
based on that of |16| and solves the Sasaki-Nakamura equation [l7| using Burlisch-Stoer integration. The accuracy 
level of the code is at the level of 10 _6 -10" 4 . 1 

First indication to the increase in the number of necessary fc-modes with increasing eccentricity was shown in [f| . 



It should be noticed, however, that the results presented in [6j for e = 0.8 may be more indication of the failure of the 
solution of the radial Teukolsky equation in the frequency domain and the need for the Sasaki-Nakamura formulation 
thereof, than a true behavior of the fc-modes. Figure [8] shows the flux in each fc mode for m = 2, after all £-modes were 
summed over, for a/M = 0.9 and p/M = 4.64 for different values of the eccentricity. Similar figures were obtained also 
for other values of m. Figures [8] and [9] suggest the following types of behavior: first, the value of fc for which the flux 
is maximal shifts to higher values with the increase in eccentricity of the orbit. Second, with increasing eccentricity, 
the flux curve broadens (mostly for positive values of fc, so that the curve becomes less and less symmetric in fc for 
high e). We also find (Fig. fTU)) that the peak of the flux curves become m-modal, i.e., at high values of e the flux 
curve breaks into a number of peaks equal to the mode m. Similar behavior was reported first in Q. The peaks 
background exhibits interesting oscillations, that appear to be independent (or, at the most, weakly dependent) of 
to, e. The broadening of the flux curves suggests that more fc-modes are required in order to obtained the total flux. 
In what follows we quantify this statement and make it precise. 



Recently, possible inaccuracies for very large k values were noticed. These inaccuracies occur for larger values of k than those used here, 
and therefore do not affect our results. 



III. SUMMATION OVER MULTIPOLES I AND k MODES 
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FIG. 4: The number of m-modes necessary for the sum over m-modes to have a certain accuracy. We plot contour curves 
at a fixed level of accuracy, showing the number of m-modes needed as a function of the eccentricity e. The contours are at 
accuracy levels of 10" 1 , 10" 2 , 10" 3 , 1CT 4 , 1CT 5 , and 1CT 6 . 

The number of fc-modes required at a fixed value of the eccentricity e to obtain the total flux to a desired accuracy 
increases logarithmically with the accuracy. This behavior is described in Fig. [11] that plots the error involved in the 
inclusion of the Nk modes (of the greatest partial fluxes) as a function of Nk for various values of the eccentricity 
e. In all cases this behavior is exponential. Notice, however, that the slope of the fitted curve becomes steeper with 
increasing eccentricity. This behavior suggests that as the eccentricity increases, the rate at which the number of 
fc-modes that are required in order to obtain the same level of accuracy increases. This rate also increases with the 
mode number m (Fig. 1 12)) . Notice, from Fig. 1 101 that also the number of negative fc- modes increases with the mode 
number m (for fixed eccentricity), although not as fast as the number of positive fc-modes. In Fig. Q2] we show the 
number of fc-modes required for a given accuracy as a function of the eccentricity e, for different values of the 
mode number m. The increase in Nk is exponential in the eccentricity e. For a given accuracy level, the exponential 
increase with e becomes steeper as the mode number m increases, as is shown in Fig. Q31 

IV. MODES COMPUTATION TIME 

The computation time of a single fc mode is a function of k. To make a prediction of the fc dependence of the 
computation time for a fc mode we first write the radial Teukolsky equation 

A 2 ± (V 1 V (r) R imu) (r) = -T tmu >(r) (1) 

where the potential V(r) is a complex valued function of r and the angular frequency to [15]. The source term 
Te mtJ (r) is a certain known function of the stress energy of the particle and its trajectory, and of the background 
geometry. As we are interested here only in the qualitative features of the solutions for Eq. (JTJ), the details of the 
source are unimportant for us here. We next focus attention on the far field limit of Eq. JT]) . To gain a qualitative 
understanding of the structure of the solution, we next consider only the asymptotic solution as r — -» oo. In that 
limit, R(r — > oo) ~ r 3 e iujr * (corresponding to an outgoing solution), where is the usual Kerr spacetime tortoise 
coordinate defined by dr*j dr = (r 2 + a 2 )/A. The solution is an oscillatory solution in r with angular frequency u>. 
The typical length scale over which the solution is oscillating is therefore A ~ 2h/lo: the greater u, the shorter the 
distance over which the solution oscillates. 

The method we use to solve the frequency-domain equation is by a Burlisch-Stoer algorithm [l8j , that successively 
divides a single step into many substeps, until a (polynomial or rational) interpolation is accurate enough. If the 
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FIG. 5: The real part of the waveform on the equatorial plane at r/M = 625 for a/M = 0.9, p/M = 4.64, and e = 0.1 for ip m =2 
(dotted), -02 + ip3 (dash-dotted), tp2 + tp3 + V>4 (dashed) and ip2 + 4>3 + ^4 + ips (solid), as a function of t/M. 
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FIG. 6: The imaginary part of the waveform on the equatorial plane at r/M — 625 for a/M = 0.9, p/M = 4.64, and e — 0.1 
for ip m =2 (dotted), ip2 + ^3 (dash-dotted), ip2 + ^3 + ip4 (dashed) and ip2 + tps + ipA + ips (solid), as a function of t/M. 



original step happens to be too large (i.e., after a certain pre-determined substep is calculated and the required 
accuracy is still not obtained), the original step will be bisected, and the process proceeds. Therefore, the greater oj 
and the shorter the oscillation length scale, the more divisions of the initial step are required, and the greater the 
number of substeps used to find the solution. As the angular frequency u> m k = mfi^ + kQ r , it is clear that the larger 
k (for a fixed value of m), the shorter the distance scale A. We therefore expect the computation time of a single k 
mode to increase with k. In fact, the Burlisch-Stoer algorithm suggests that as Lo m k increases, the more substeps are 
needed, until for some value of k dividing the step into substeps is no longer sufficient, and the step is bisected. Then, 
division into substeps is again sufficient, until another k value is obtained for which the step is bisected. Each time 
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FIG. 7: The phase of the waveform on the equatorial plane for the same parameters as in Fig. [5] Shown are the phases 
corresponding to ip m =2 (x), tJj 2 + 4>3 (□), 4>2 + ip3 + ip4 (*) and ip2 + ip3 + ip4 + ips (°), as a function of t/M. 



the step is bisected, the total number of substeps computed jumps discontinuously, so that we expect the k mode 
computation time to increase as a staircase function. In Fig. [15] we show the computation time of an individual k 
mode as a function of k, and also the number of Burlisch-Stoer iterations done. They both behave as expected. The 
preceding discussion suggests that a similar behavior of the computation time is found also as a function of m. This 
is indeed seen in Fig. [TH 

As can be seen from Figs. [TS] and [TH] the discontinuous jumps in the mode computation times occur at different 
values of m, k. Indeed, our discussion about explains this behavior: jumps occur when the angular frequencies uj m k 
arrive at some threshold values, corresponding to threshold wavelengths. As these threshold values should be about 
the same, we expect jumps to occur for mj,kj satisfying kj£l r ~ m,j$l v . Notice that kj (rrij) is defined here for 
vanishing m (k). For the data presented here, for e = 0.4 we find kj/rrij w 1.33, while H, ip /il r w 1.52; for e = 0.5 
we find kj/rrij « 1.40, while Vl v /VL r k, 1.47; and for e — 0.6 we find kj/rrij « 1.44, while Vt v /Q, r ps 1.39, so that the 
difference between kj/rrij and fl v /Q r is at order 10%. We argue that these results support our interpretation of the 
discontinuous jumps in the computation time. 

Our discussion suggests that for a fixed value of k (m) and varying m (fc), there are values of the latter for which 
the computational time drops. These are the values that approximately satisfy Lu m k ^ £l v , fi r - (Recall that fc, m can 
be either positive or negative.) That is, while fixing either k or m and varying the other, we expect to find a drop in 
the computation time (or, equivalently, in the number of Burlisch-Stoer iterations.) Indeed, we find this behavior, as 
is shown in Fig. [T7] Notice, that for k = the drop in computational time is found for m = 0, and as k increases, so 
does the value of m for which the drop is found. We also comment that this drop is broadened with increasing values 
of k. 

The dependence of the computational time on the mode number k is therefore rather intricate. It can be approxi- 
mated as follows: For each value of the eccentricity e we take in practice the computational time of a single mode k 
to be a linear function of k, that fits the computational data. This approximation underestimates the computation 
time for some k modes, specifically those immediately following a discontinuous jump in the computation time, and 
overestimates the computation time for k modes just before a jump. However, as we are mostly interested in the total 
computation time of the sum over all modes, this approximation may be quite reasonable for the sum over all modes. 
In this approximation the sum over k modes up to some fc max is a quadratic function of fc max . Notice, that we can 
use Fig. llll and ll2l to find how many k modes we need to sum over to obtain the required accuracy level. 
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FIG. 8: The flux in different fc-modes, after summation over all £, for m — 2, for a/M = 0.9 and p/M = 4.64: e = 0.1 (dotted), 
e = 0.2 (dash-dotted), e = 0.3 (thin dashed), e = 0.4 (thin solid), e = 0.5 (thick dashed), and e = 0.6 (thick solid). For each 
value of the eccentricity, the flux is normalized to the flux at the value of k for which the flux is maximal. 
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FIG. 9: Same as Fig. HI for m = 3. 

V. COMPARISON OF TOTAL COMPUTATION TIME OF TD AND FD CALCULATIONS 

Finally, we put together all the previous results, to compare the actual total calculation time of TD and FD codes. 
We choose an orbit that gives no special preference for one approach over the other, i.e., a moderate value for the 
eccentricity. In practice, we take e = 0.5 and p = 4.6AM. 

For the FD calculation we first set the desired level of accuracy, which in practice we choose here to be 10 -3 , and 
then use Fig. [4] (or the data in Table Q} to determine the number of m modes one needs to sum over to guarantee the 
desired accuracy level. The needed number of m modes is found in practice to be 7, for m = 2, 3, • • • , 9. 

For each individual m mode we then determine the number of k modes necessary to obtain the desired accuracy 
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FIG. 10: Same as Fig. [8] for m = 2 (dotted), m — 3 (dashed), and m = 4 (solid), for e = 0.6. Notice that the first two curves 
repeat data already presented in Figs. [8] and [9] 



level. This determination is done with data as in Figs. [I3l and [Til The values of the k modes we used in practice are 
as follows: for m = 2 we used —2 $C k ^ 10, for m = 3 we used —3 ^ k ^ 13, for m = 4 we used —3 ^ k ^ 16, for 
rn = 5 we used — 2 ^ fc ^ 21, for m = 6 we used — 1 ^ k ^ 23, for to = 7 we used — 2 ^ fc ^ 25, for m = 8 we used 
— 1 ^ k ^ 28, and for m = 9 we used — 1 ^ k ^ 33. In all cases we computed 2 ^ £ ^ 5. The data are summarized 
in Table |TT] and in Fig. [T9j The FD calculations were optimized for most efficient calculation for the desired accuracy 
level. Improvements to the computation time can be achieved, but only at the level of the code, e.g., making the 
Burlisch-Stoer algorithm more efficient. We believe that while such improvements can be made, their effect would be 
moderate. In this sense, the FD curve in Fig. [TH] represents a lower bound on the FD calculation time. Increasing the 
eccentricity of the orbit e would result in the FD curve of Fig. [l"9l moving up in the figure (increase in the calculation 
time of each m mode, because of increasing number of required number of k modes and increase in the computation 
time of individual k modes) and also having a faster growth rate. 



TABLE II: Number of k modes for each m modes in the FD calculation, and the corresponding computation time in both the 
FD and TD calculations. The orbital parameters are e = 0.5 and p/M = 4.64 and a/M — 0.9. 



m mode 


Range of k modes 


Total k modes 


FD time (sec) 


TD time (sec) 


2 


-2 k 10 


13 


51 


954 


3 


-3 s; k < 13 


17 


99 


945 


4 


-3 < k < 16 


20 


116 


952 


5 


-2 < k < 21 


24 


179 


942 


6 


-1 < k < 23 


25 


218 


953 


7 


-2 < k < 25 


28 


254 


953 


8 


-1 ^ k < 28 


30 


303 


952 


9 


-1 ^ k < 33 


35 


391 


942 



For the TD calculations we used the 2+ ID code of [j], with radial resolution of Ar = M/20 (and temporal 
resolution of At = Ar/2), and angular resolution of A9 = tt/32. We placed the inner and outer boundaries at 
r*/M = —50 and at r*/M — 350, respectively. We approximate the scattering problem on the boundaries by 
assuming Sommerfeld boundary conditions, i.e., no incoming radiation from outside the boundaries. One may of 
course make this approximation better by pushing the outer boundary outwards (and the inner boundary inwards). 
Reflections from the boundaries then do not contaminate the center of the computational domain (at r*/M = 150) 
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FIG. 11: The error A in summing over a number Nk of fc-modes, after summation over all I, for m = 2, for a/M — 0.9 and 
p/M = 4.64: e = 0.1 (o) , e = 0.2 (*), e = 0.3 (x), e = 0.4 (□), e = 0.5 (o), and e = 0.6 (A). For each value of the eccentricity, the 
solid line describes a fitted line as follows, corresponding to increasing eccentricity: log 10 A = 1.2796- 1.1630 N k {R 2 = 0.9810), 
log 10 A = 2.2131 - 0.9567 N k (R 2 = 0.9691), log 10 A = 2.5781 - 0.7740 N k (R 2 = 0.9774), log 10 A = 3.0085 - 0.6445 iV fe 
(R 2 = 0.9962), log 10 A = 3.2273 - 0.5197 N k {R 2 = 0.9987), log 10 A = 3.8898 - 0.4267 N k (R 2 = 0.9383). For each fitted curve 
we include the square of the correlation coefficient R 2 




FIG. 12: The error A in summing over a number N k of fc-modes, after summation over all £, for a/M = 0.9, p/M — 4.64, 
and e = 0.3: m = 2 (o), m = 3 (□), and m — 4 (*). For each value of the eccentricity, the solid line describes a fitted line 
as follows, corresponding to increasing value of m: log 10 A = 2.5781 - 0.7740 N k (R 2 = 0.9774), log 10 A = 2.7038 - 0.6366 N k 
(R 2 = 0.9950), log 10 A = 3.9493 - 0.6126 N k (R 2 = 0.9893). For each fitted curve we include the square of the correlation 
coefficient R 2 
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FIG. 13: The number N k of fc-modes required for a given accuracy as a function of the eccentricity e. Shown are 6 levels 
of accuracy, and their corresponding contours: 1CT 1 (o), 1CT 2 (*), 1CT 3 (x), 1CT 4 (□), 1CF 5 (o), and 10~ 6 (A), for three 
values of the mode number m = 2,3,4. Upper panel (m = 2): the solid lines describe a fitted line as follows, corresponding 
the increasing accuracy: log 10 N k = 0.2253 + 1.3355 e (R 2 = 0.9642), log 10 N k = 0.3632 + 1.2783 e (R 2 = 0.9872), log 10 N k = 
0.4671 + 1.2442 e (R 2 = 0.9945), log 10 N k = 0.5506 + 1.2215 e {R 2 = 0.9959), log 10 N k = 0.6206 + 1.2053 e (R 2 = 0.9948), 
log 10 N k = 0.6808+1.1930 e (R 2 — 0.9928). Middle panel (m = 3): the solid lines describe a fitted line as follows, corresponding 
the increasing accuracy: log 10 N k = 0.3040 + 1.1118 e (R 2 = 0.8860), log 10 N k = 0.4321 + 1.2321 e (R 2 = 0.9727), log 10 N k = 
0.5329 + 1.2886 e (R 2 = 0.9900), log 10 N k = 0.6151 + 1.3216 e (R 2 = 0.9934), log 10 N k = 0.6844 + 1.3433 e (R 2 = 0.9932), 
log 10 N k = 0.7442 + 1.1588 e (i? 2 = 0.9919). Lower panel (m = 4): the solid lines describe a fitted line as follows, corresponding 
the increasing accuracy: log 10 N k = 0.4636 + 0.8385 e (R 2 = 0.7739), log 10 N k = 0.5159 + 1.2410 e (i? 2 = 0.9747), log 10 A^ fe = 
0.5841 + 1.4111 e (i? 2 = 0.9894), log 10 N k = 0.6478 + 1.5066 e (i? 2 = 0.9884), log 10 JV fc = 0.7050 + 1.5681 e (i? 2 = 0.9853), 
log 10 N k = 0.7564 + 1.6111 e (i? 2 = 0.9821). 

40 | 1 1 1 1 1 1 1 1 1 1 1 1 




1 1 1 1 1 1 1 1 1 1 1 1 

0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 

e 

FIG. 14: The number iVfc of fc-modes required for an accuracy of 10 -4 as a function of the eccentricity e. Shown are 3 values 
of m, and their corresponding contours, based on a fit to an exponential: m = 2 (o), m = 3 (*), and m = 4 (x). 
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FIG. 15: The dependence of the computation time of a k mode on k: the computation time t (normalized by the maximal 
computation time is shown in x , and the corresponding normalized number of iterations n,: that the Burlisch-Stoer engine does 
is shown in o. The data are taken for m = 0, and for p/M = 4.64 and e — 0.4. 
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FIG. 16: Same as Fig. [15] for the m modes. Here, k — and £ = 20. 



until t/M = 400, and we integrate in time until then. In practice, we wait at each evaluation point until the initial 
spurious waves (that result from imprecise initial data that correspond to a particle suddenly appearing at t/M = 0), 
and then average the flux over a full period of the orbit. The angular period is = 2n Af -1 / 2 (r 3 / 2 +aM 1 ^ 2 ), where 
r is the semimajor axis. For our choice of parameters, w 102. AM. The radial period is found to be T r rs 162. 5M, 
and it is the latter (radial) period we use for our analysis. 

We extract the (radial period averaged) flux at a number of extraction distances, in practice at r*/M = 
30,40, • • • , 100, and then fit the (finite extraction distance) fluxes to an inverse-square function as in Eq. (|A2j) . This 
allows us to achieve fluxes that agree with the FD fluxes to 10~ 3 , while not having to integrate to very late times (and 
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FIG. 17: Number of iterations of the Burlisch-Stoer solver as a function of m for fixed values of k, for the same orbital 
parameters as above, for e = 0.1. 




0.1 



FIG. 18: The total computation time of the sum over k modes from k — to k = fc max as a function of the eccentricity e. The 
computation time is presented in units of the maximal computation time in the chosen range of parameters. 



commensurately increase also the spatial computational domain). Notably, we can obtain higher accuracy than that 
reported in [tJ while integrating to a shorter time because we make use of the finite time corrections to the energy 
flux. We discuss this method in detail in the Appendix. 

The TD computational time is presented in Table ITT1 It can be made more efficient using code improvements, and 
data analysis improvements. First, we used gaussian source modeling as in [7] , and the discrete 5 approach may 
improve efficiency considerably. In fact, Ref. 0] has argued for a full order of magnitude reduction in computation 
time. While this statement may be an overestimate of the numerical capabilities — especially for generic (i.e., eccentric 
or inclined) orbits — it is certainly possible to improve the efficiency of the TD code by changing the method of 
calculation of the source term. (The discrete delta approach has also other advantages over the gaussian model. See 
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8] for more detail.) An important property of the TD code is that it is very naturally parallelizable. Certainly the FD 
computation is also parallelizable, most naturally by having each k mode computed on a different processor. In such a 
case, the total FD computation time is controlled by the k mode that takes the longest to compute, i.e., the largest k 
(see Fig. 1151) . Another improvement of efficiency of the FD code may come from improvements of the Burlisch-Stoer 
engine, specifically the solution for logr instead of r. It is currently unknown how the computational efficiency and 
numerical accuracy would be affected by such a change. 

One factor that limits a possible reduction of the TD computation time is the need to average over an orbital 
(radial) period to find the average flux per period. Even for the strong field orbit considered here, with p/M = 4.64 
and e = 0.5, the radial period is T r = 162. 5M, which puts a lower bound on how short the computational integration 
can be. One could perhaps shorten the computation time if only half an orbital period is computed, but because the 
phase of the orbit at the end of the spurious wave epoch is arbitrary, on the average at least three quarters of an orbit 
need to be computed. However, we were interested here in comparison of the fluxes to infinity. If one is interested 
in the wave function itself, shorter computation times can be used, thus saving much of the TD computation time. 
No equivalent save in time can be achieved with the FD computation. Reduction of the TD computation time would 
result in the dashed line in Fig. [19] moving lower, thus making the computational time of TD and FD comparable 
already for lower m values. Most importantly, when higher eccentricity orbits are considered, the TD computation 
time remains unchanged, whereas the FD computation time increases considerably. More specifically, not only is the 
solid curve in Fig. rj5] shifted up, it also becomes steeper (i.e., its rate of growth increases), so that the TD and FD 
computation times become comparable for lower m values. In addition, the TD calculation does not require any extra 
computation to find the waveform, as it is already available. In fact, an extra computation was needed to find the 
fluxes. On the other hand, the FD calculation requires an extra computation to produce the waveform, which we 
have not done here. Therefore, the greater efficiency of FD over TD for the modes shown here is much less impressive 
when the computation time for waveforms is of interest. 

We therefore suggest that at high eccentricity orbits the computational efficiency of the TD and the FD approaches 
is comparable for high m values. We believe that the question of the total (sum over m) computation time is not 
necessarily a very important one, and suggest that for actual computations the low m modes are calculated using the 
FD approach, and the high m modes are calculated using the TD approach. The determination of which m values are 
low and which arc high depends of course on the parameters of the system, and the accuracy level of the computation. 
Using a single method, however, could still possibly make TD more efficient than FD for very high values of the 
eccentricity. 

In the analysis of this paper we have focused on equatorial orbits. Such orbit simplify the problem not just 
computationally: the number of dimensions in the FD parameter space is significantly lower than that for generic 
orbits. While the analysis of equatorial orbits presented here is non-trivial, the added complications of handling generic 
orbits warrant separate treatment. When generic orbits are considered, finding the number of i modes required to 
achieve a set accuracy level with the FD approach becomes a non-trivial problem, as is the question of finding the 
computational time of each £ mode as a function of the mode number. Notably, the TD calculations remain unchanged, 
as is the associated computation time. We therefore expect the FD computation time to increase significantly with 
the transition from equatorial to inclined and then generic orbits, an increase that has no counterpart in the TD 
case. To do such an analysis, one needs to check the behavior of all combinations of the modes fc, £, m, find the 
minimal number of modes required for a set accuracy level, and then compute those modes. The TD calculation 
requires that the same number of m modes are computed, and this computation time depends only on the duration 
of the computation, the size of the computational domain, and of course the grid resolution. Modest improvements 
to the TD computation time may be achieved by improving the efficiency of the source calculations (e.g., the discrete 
delta approach), but they are not expected to be as dramatic as previously suggested Q. We expect the threshold 
eccentricity above which the TD approach is more computationally efficient than the FD approach for some of the 
needed m modes to be moderate-high, specifically in the range of e ~ 0.6-0.7. 
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FIG. 19: Comparison of computation time for FD and TD for an orbit with p/M = 4.64 and e = 0.5 as a function of the mode 
number m. The circles (o) denote FD computation time, and the squares (□) show TD computation times. Both sets of runs 
were done on the same computer, in this case a single processor of a dual 2.7GHz PowerPC G5. The red curve is the fitted 
curve t — 21.471 m 1,2919 sec, that fits the numerical data with a squared correlation coefficient of R 2 = 0.9895. For both cases 
the accuracy level is set at 10 -3 . 



APPENDIX A: NOTES ON WAVE-EXTRACTION DISTANCE FOR THE TIME DOMAIN CODE 



The peeling theorem describes how the various components of the gravitational field of a bound gravitational system 
behave as the observer moves away from the source (see, e.g., fl9l . |20|). Specifically, in the wave zone of a radiating 
system, the "J 4 component of the Weyl curvature — which describes in the wave zone the outgoing transverse radiative 
degrees of freedom — typically drops off as 'J 4 ~ r _1 , where r is the distance from the radiating system. The peeling 
theorem also gives the drop off rates for the remaining Weyl scalars, and also for spin coefficients and other Newman- 
Penrose quantities. In addition to the leading order decay rate, the peeling theorem has been extended also to include 
correction terms in r _1 . Specifically, Newman and Unti showed that [2lT | 



* 4 = — - 



^0 2a°^l 



OfcVTrO 

T 3 ,k 



+ 0(r- J ), 



where a is a spin coefficient and £ z is an "integration constant" from the radial equations. A superscript "0" means 



the coefficient of the slowest decaying term in an expansion in 



As the flux of energy in gravitational waves at 



infinity T ~ lim^oo r 2 ^\ ~ (^T) 2 > ^4 is directly related to the flux at infinity. 

Notably, the leading order corrections in r -1 to ^4 are proportional to the asymptotic value of the Weyl scalar ^3 
and its gradient. In the Teukolsky formalism, that conveniently describes black hole perturbations, the Weyl scalar 
^3 = ("transverse frame"), so that the leading order correction to ^4 is at the next order: ^4 — ^4/r + 0(r -3 ). 
Therefore, 



T ~ lim 



r 2 * 2 ~ lim r 2 

r — >oo 



+ 0(r~ 3 ) 



(*°) 2 [l + 0(r- 2 )]. 



(Al) 



We therefore are motivated to introduce the ansatz 

E = E 00 [l-A(mX/r) 2 



(A2) 



where A : 
mX = f2" 



A/(2tt) = (rl /2 + aM l / 2 )/(mM x l 2 ), and r is the Boyer-Lindquist radius of the orbit. Notice that 
so that our ansatz can be written as E = i?oo[l — A(flr)~ 2 ]. Note that the expansion parameter 
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(fir) 2 <C 1, as the field is evaluated far from the radiating system, specifically outside the "light cylinder." This 
ansatz is also suggested by Fig. 1 of Ref . jjj and by Table VI of Ref. [H , where a best fit was done to the ansatz 
-Eskh = -Eoo[l — <l{ r o/ r ) v ]i and the free parameter p was found numerically to be rather close to 2 (with deviations of 
< 3%) for a large range of orbits. 

We test our ansatz by fitting the outgoing flux of energy for a number of circular and equatorial orbits around a 
Kerr black hole as detected at a sequence of distances from the radiating system to the ansatz (|A2|I . In Table IIIII we 
show the outgoing fluxes from particles in circular orbits around a Kerr black hole, taken at a sequence of distances. 
Table IIVI shows the squared correlation coefficient for the various cases we have checked. In all cases the high value 
of the correlation coefficient corroborates our ansatz. In Table fVl we show the values of the free parameter A for the 
same cases. Lastly, in Table IVT1 we compare the flux to infinity based on the ansatz (|A2[) to the flux obtained from a 
frequency-domain calculation. Notably, the associated errors appear to be larger for m = 4, 5 than for m = 2,3. We 
attribute these higher errors to the lower values of the flux for high values of m. Specifically, to obtain the flux with 
the time domain one needs to first subtract the "flux" due to the spurious radiation associated with the initial time 
of the simulation Q. As the fluxes become smaller, this subtraction becomes less accurate. 

TABLE III: Fluxes per unit mass extracted from the time-domain code at a sequence of radii on the numerical grid. En is the 
flux per unit mass measured at a radius RM. For |m| =2,3 the data are identical to table V of Ref. Notice that in IsH 
(and also tables 1,11 of 0]) each value of \m\ is in fact the sum of the contributions of m and — m. 



\m\ 


r /M 


a/M 




E200 


-E300 


-B400 


-B500 


£-600 


1 


4.0 


0.99 


1.2779 x 10~ 6 


1.3049 x 10~ 6 


1.3107 x 10~ 6 


1.3130 x 10~ 6 


1.3143 x 10~ 6 


1.3150 x 10~ 6 


2 


4.0 


0.99 


1.2284 x 10" 3 


1.2341 x 10" 3 


1.2351 x 10" 3 


1.2355 x 10~ 3 


1.2356 x 10~ 3 


1.2357 x 10" 3 


3 


4.0 


0.99 


2.9481 x 10~ 4 


2.9639 x 10~ 4 


2.9667 x 10~ 4 


2.9677 x 10~ 4 


2.9681 x 10" 4 


2.9682 x 10~ 4 


4 


4.0 


0.99 


8.3615 x 10" 5 


8.4525 x 10 -5 


8.4769 x 10 -5 


8.4865 x icr 6 


8.4924 x 10~ 5 


8.4955 x 10" 5 


5 


4.0 


0.99 


2.6586 x 10~ 5 


2.6875 x 10~ 5 


2.6949 x 10~ 5 


2.6982 x 10~ 5 


2.7000 x 10~ 5 


2.7015 x 10~ 5 


1 


10 


0.90 


1.9565 x 10" 8 


2.4481 x 10" 8 


2.5424 x 10~ 8 


2.5735 x 10~ 8 


2.5888 x 10~ 8 


2.5968 x 10~ 8 


2 


10 


0.90 


2.0865 x 10~ 5 


2.1965 x 10~ 5 


2.2161 x 10 -5 


2.2228 x 10~ 5 


2.2259 x 10~ 5 


2.2275 x 10~ 5 


3 


10 


0.90 


2.3396 x 10" 6 


2.4794 x 10" 6 


2.5043 x 10" 6 


2.5128 x 10 -6 


2.5167 x 10~ 6 


2.5188 x 10 -6 


4 


10 


0.90 


3.2046 x 10~ 7 


3.4170 x 10~ 7 


3.4576 x 10~ 7 


3.4723 x 10 -7 


3.4796 x 10~ 7 


3.4836 x 10~ 7 


5 


10 


0.90 


4.8313 x 10" 8 


5.1497 x 10" 8 


5.2112 x 10" 8 


5.2330 x 10~ 8 


5.2438 x 10" 8 


5.2500 x 10 -8 


1 


10 


0.99 


1.6454 x 10 -8 


2.0725 x 10~ 8 


2.1513 x 10 -8 


2.1778 x 10~ 8 


2.1907 x 10~ 8 


2.1976 x 10~ 8 


2 


10 


0.99 


2.0516 x 10~ 5 


2.1605 x 10" 5 


2.1799 x 10" 5 


2.1884 x 10 -5 


2.1914 x 10~ 5 


2.1931 x 10~ 5 


3 


10 


0.99 


2.2889 x 10~ 6 


2.4279 x 10~ 6 


2.4526 x 10~ 6 


2.4610 x 10~ 6 


2.4629 x 10~ 6 


2.4670 x 10~ 6 


4 


10 


0.99 


3.1481 x 10" 7 


3.3596 x 10" 7 


3.3998 x 10" 7 


3.4144 x 10 -7 


3.4208 x 10~ 7 


3.4329 x 10~ 7 


5 


10 


0.99 


4.7238 x 10~ 8 


5.0395 x 10" 8 


5.1012 x 10~ 8 


5.1174 x 10~ 8 


5.1122 x 10~ 8 


5.1182 x 10~ 8 


1 


12 


0.0 


1.9126 x 10~ 8 


2.8001 x 10~ 8 


2.9932 x 10~ 8 


3.0499 x 10~ 8 


3.0772 x 10~ 8 


3.0928 x 10~ 8 


2 


12 


0.0 


0.9792 x 10~ 5 


1.0628 x 10~ 5 


1.0777 x 10~ 5 


1.0827 x 10~ 5 


1.0850 x 10~ 5 


1.0862 x 10~ 5 


3 


12 


0.0 


0.9769 x 10" 6 


1.0663 x 10" 6 


1.0825 x 10" 6 


1.0881 x 10~ 6 


1.0906 x 10 -6 


1.0920 x 10~ 6 


4 


12 


0.0 


1.1883 x 10 -7 


1.3096 x 10" 7 


1.3321 x 10~ 7 


1.3401 x 10~ 7 


1.3440 x 10~ 7 


1.3462 x 10~ 7 


5 


12 


0.0 


1.5953 x 10" 8 


1.7590 x 10" 8 


1.7874 x 10" 8 


1.7982 x 10 -8 


1.8033 x 10~ 8 


1.8062 x 10 -8 



TABLE IV: Correlation coefficient R 2 for the Ansatz. The roman numerals list the four orbits considered. I: ro = AM, 
a = 0.99M. II: r = 10M, a = 0.9M. Ill: r = 10M, a = 0.99M. IV: r = 12A/, a = 0. Notice we have changed the order of 
the orbits compared with previous paper, to arrange them in ascending order of wavelengths. For all cases, the fit is done with 
6 data points, at values of r/M — 100, 200, 300, 400, 500, 600. (The case 11,4 has a different extraction locations.) 



m mode 


1 


2 


3 


4 


5 


I 


0.998976026 


0.99998427 


0.999899201 


0.9956862 


0.99507507 


II 


0.99999139 


0.99996189 


0.99996031 


0.999946517 


0.99994224 


III 


0.99999644 


0.99996171 


0.99995783 


0.998752311 


0.99997298 


IV 


0.99978769 


0.99994700 


0.99997873 


0.99999386 


0.99997505 
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TABLE V: Best-fit value for the parameter A. The roman numerals designate the same orbits as in Table ILTTl The bracketed 
number denotes the uncertainly in the last figure. 



m mode 


1 


2 


3 


4 


5 


I 


3.5458(2) 


0.75753(2) 


0.86647(5) 


1.9655(3) 


1.9622(4) 


II 


2.3812(2) 


0.6156(4) 


0.6915(4) 


0.7751(5) 


0.7716(5) 


III 


2.4132(1) 


0.6224(4) 


0.6976(4) 


0.799(2) 


0.7734(5) 


IV 


2.250(2) 


0.5857(7) 


0.6287(5) 


0.6950(3) 


0.6917(5) 



TABLE VI: Comparison of frequency-domain fluxes with time-domain fluxes extrapolated to infinty based on the ansatz ()A2|| . 



\m\ 


r /M 


a/M 


Eoo 


-Bfd 


{Eoo — Efd) 1 Efd 


1 


4.0 


0.99 


1.3153 x 10~ 6 


1.3403 x 10~ 6 


-0.0187 


2 


4.0 


0.99 


1.2359 x 10" 3 


1.2418 x 10" 3 


-0.0047 


3 


4.0 


0.99 


2.9689 x 10~ 4 


2.9621 x 10~ 4 


0.0023 


4 


4.0 


0.99 


8.4995 x 10" 5 


8.6330 x 10" 5 


-0.0155 


5 


4.0 


0.99 


2.7024 x 10~ 5 


2.7162 x 10~ 5 


-0.0051 


1 


10 


0.90 


2.6147 x 10~ 8 


2.6555 x 10" 8 


-0.0154 


2 


10 


0.90 


2.2320 x 10" 5 


2.2281 x 10" 5 


0.0018 


3 


10 


0.90 


2.5245 x 10~ 6 


2.5221 x 10" 6 


0.0010 


4 


10 


0.90 


3.4902 x 10 -7 


3.5345 x 10" 7 


-0.0125 


5 


10 


0.90 


5.2599 x 10~ 8 


5.3255 x 10~ 8 


-0.0123 


1 


10 


0.99 


2.2138 x 10" 8 


2.2503 x 10" 8 


-0.0162 


2 


10 


0.99 


2.1969 x 10~ 5 


2.1974 x 10~ 5 


-0.0002 


3 


10 


0.99 


2.4726 x 10" 6 


2.4709 x 10" 6 


0.0007 


4 


10 


0.99 


3.4388 x 10~ 7 


3.4467 x 10~ 7 


-0.0023 


5 


10 


0.99 


5.1469 x 10" 8 


5.1687 x 10" 8 


-0.0042 


1 


12 


0.0 


3.1227 x 10~ 8 


3.1456 x 10~ 8 


-0.0073 


2 


12 


0.0 


1.0897 x 10" 5 


1.0861 x 10" 5 


0.0033 


3 


12 


0.0 


1.0956 x 10 -6 


1.0945 x 10~ 6 


0.0010 


4 


12 


0.0 


1.3503 x 10" 7 


1.3658 x 10" 7 


-0.0114 


5 


12 


0.0 


1.8121 x 10~ 8 


1.8317 x 10 -8 


-0.0107 
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